### CREATED BY JESSICA SCHOENEHRR AND NICK WATERBURY ###
### Replication files for "Confessions at the Supreme Court" ###
### File 3 - Replication of Appendix Models and Figures ###

### data sources:
### - Confessions of error originally collected by Bruhl (2009)
### 		- Updated by authors through the 2014 term
### - Politicization measure from Patrick Wohlfarth (2009)
###		- Updated by the authors through the 2014 term
### - sgJusticeDistance from the Judicial Common Space (Epstein, Martin, Segal, and Westerland)
###		- http://epstein.wustl.edu/research/JCS.html
### - CSI from Collins and Cooper
###		- https://dataverse.harvard.edu/dataset.xhtml?persistentId=doi:10.7910/DVN/UR2KYE
### - SG list of cases and SG as petitioner from Black and Owens (2012)
###		- updated by the authors per the OSG website
###		- https://www.justice.gov/osg/supreme-court-briefs
### - all other data from the Supreme Court DataBase (Spaeth, Epstein, et al.)
### 		- http://scdb.wustl.edu/index.php

library(arm)
library(compactr)
library(stargazer)
library(haven)

library(miceadds)
library(faraway)

library(tidyverse)
library(ggplot2)
library(ggthemes)
library(gridExtra)
library(cowplot)
library(modelsummary)

data1 <- read_dta("MeritsJusticeCentered20201026.dta")

data1$justice <- as.factor(data1$justice)
data1$votesWithSG <- as.factor(data1$votesWithSG)

#############
### MODEL ###
#############

model4MeritsFullMLM <- glm.cluster(votesWithSG ~ policyChange
							+ noPolicyChange
							+ coeNoDetail
							+ coeCriminalCase
							+ coePastTally
							+ coeInPrevious
							+ politicization
							+ sgJusticeDist
							+ sgPetitioner
							+ constitutionalCase
							+ declareUncon
							+ CSI
							+ realOutlierSignal
							+ civilRightsCase,
							data = data1,
							cluster = "justice",
							family = binomial(link = logit))
							
summary(model4MeritsFullMLM$glm_res)
logLik(model4MeritsFullMLM$glm_res)
BIC(model4MeritsFullMLM$glm_res)
nobs(model4MeritsFullMLM$glm_res)

data2 <- read_dta("CVSGJusticeCentered20201026.dta")

data2$justice <- as.factor(data2$justice)
data2$votesWithSG <- as.factor(data2$votesWithSG)

#############
### MODEL ###
#############

model4CVSGFullMLM <- glm.cluster(votesWithSG ~ policyChange
							+ noPolicyChange
							+ coeCriminalCase
							+ coePastTally
							+ coeInPrevious
							+ politicization
							+ sgJusticeDist
							+ sgPetitioner
							+ constitutionalCase
							+ declareUncon
							+ CSI
							+ realOutlierSignal
							+ civilRightsCase,
							data = data2,
							cluster = "justice",
							family = binomial(link = logit))
							
summary(model4CVSGFullMLM$glm_res)
logLik(model4CVSGFullMLM$glm_res)
BIC(model4CVSGFullMLM$glm_res)
nobs(model4CVSGFullMLM$glm_res)

data3 <- read_dta("VoluntaryACBJusticeCentered20201026.dta")

data3$justice <- as.factor(data3$justice)
data3$votesWithSG <- as.factor(data3$votesWithSG)

#############
### MODEL ###
#############

model4VACBFullMLM <- glm.cluster(votesWithSG ~ policyChange
							+ noPolicyChange
							+ coeNoDetail
							+ coeCriminalCase
							+ coePastTally
							+ coeInPrevious
							+ politicization
							+ sgJusticeDist
							+ sgPetitioner
							+ constitutionalCase
							+ declareUncon
							+ CSI
							+ realOutlierSignal
							+ civilRightsCase,
							data = data3,
							cluster = "justice",
							family = binomial(link = logit))
							
summary(model4VACBFullMLM$glm_res)
logLik(model4VACBFullMLM$glm_res)
BIC(model4VACBFullMLM$glm_res)
nobs(model4VACBFullMLM$glm_res)

models1 <- list(model4MeritsFullMLM$glm_res, model4CVSGFullMLM$glm_res, model4VACBFullMLM$glm_res)
modelsummary(models1, output = "latex", stars = TRUE)

##############

###
# Drop out the pre-1992 term, include variable for coeOnDiffAdmin
###

data1Short <- data1 %>% filter(term >= 1992)

model2MeritsFullPooled <- glm.cluster(votesWithSG ~ policyChange
							+ noPolicyChange
							+ coeNoDetail
							+ coeCriminalCase
							+ coeInPrevious
							+ coePastTally
							+ politicization
							+ realOutlierSignal
							+ sgJusticeDist
							+ sgPetitioner
							+ constitutionalCase
							+ civilRightsCase
							+ declareUncon
							+ CSI
							+ civilRightsCase,
							data = data1Short,
							cluster = "justice",
							family = binomial(link = logit))
							
summary(model2MeritsFullPooled$glm_res)
logLik(model2MeritsFullPooled$glm_res)
BIC(model2MeritsFullPooled$glm_res)
nobs(model2MeritsFullPooled$glm_res)

model3MeritsFullPooled <- glm.cluster(votesWithSG ~ policyChange
							+ noPolicyChange
							+ coeNoDetail
							+ coeCriminalCase
							+ coeInPrevious
							+ coePastTally
							+ realOutlierSignal
							+ coeOnDiffAdmin
							+ sgJusticeDist
							+ sgPetitioner
							+ constitutionalCase
							+ civilRightsCase
							+ declareUncon
							+ CSI,
							data = data1Short,
							cluster = "justice",
							family = binomial(link = logit))
							
summary(model3MeritsFullPooled$glm_res)
logLik(model3MeritsFullPooled$glm_res)
BIC(model3MeritsFullPooled$glm_res)
nobs(model3MeritsFullPooled$glm_res)

###
# Drop out the pre-1992 term, include variable for coeOnDiffAdmin
###

data2Short <- data2 %>% filter(term >= 1992)

model2CVSGFullPooled <- glm.cluster(votesWithSG ~ policyChange
							+ noPolicyChange
							#+ coeNoDetail
							+ coeCriminalCase
							+ coeInPrevious
							+ coePastTally
							+ politicization
							+ realOutlierSignal
							+ sgJusticeDist
							+ sgPetitioner
							+ constitutionalCase
							+ civilRightsCase
							+ declareUncon
							+ CSI,
							data = data2Short,
							cluster = "justice",
							family = binomial(link = logit))
							
summary(model2CVSGFullPooled$glm_res)
logLik(model2CVSGFullPooled$glm_res)
BIC(model2CVSGFullPooled$glm_res)
nobs(model2CVSGFullPooled$glm_res)

model3CVSGFullPooled <- glm.cluster(votesWithSG ~ policyChange
							+ noPolicyChange
							#+ coeNoDetail
							+ coeCriminalCase
							+ coeInPrevious
							+ coePastTally
							+ realOutlierSignal
							+ coeOnDiffAdmin
							+ sgJusticeDist
							+ sgPetitioner
							+ constitutionalCase
							+ civilRightsCase
							+ declareUncon
							+ CSI,
							data = data2Short,
							cluster = "justice",
							family = binomial(link = logit))
							
summary(model3CVSGFullPooled$glm_res)
logLik(model3CVSGFullPooled$glm_res)
BIC(model3CVSGFullPooled$glm_res)
nobs(model3CVSGFullPooled$glm_res)

###
# Drop out the pre-1992 term, include variable for coeOnDiffAdmin
###

data3Short <- data3 %>% filter(term >= 1992)

model2VACBFullPooled <- glm.cluster(votesWithSG ~ policyChange
							+ noPolicyChange
							+ coeNoDetail
							+ coeCriminalCase
							+ coeInPrevious
							+ coePastTally
							+ politicization
							+ realOutlierSignal
							+ sgJusticeDist
							+ sgPetitioner
							+ constitutionalCase
							+ civilRightsCase
							+ declareUncon
							+ CSI,
							data = data3Short,
							cluster = "justice",
							family = binomial(link = logit))
							
summary(model2VACBFullPooled$glm_res)
logLik(model2VACBFullPooled$glm_res)
BIC(model2VACBFullPooled$glm_res)
nobs(model2VACBFullPooled$glm_res)

model3VACBFullPooled <- glm.cluster(votesWithSG ~ policyChange
							+ noPolicyChange
							+ coeNoDetail
							+ coeCriminalCase
							+ coeInPrevious
							+ coePastTally
							+ realOutlierSignal
							+ coeOnDiffAdmin
							+ sgJusticeDist
							+ sgPetitioner
							+ constitutionalCase
							+ civilRightsCase
							+ declareUncon
							+ CSI,
							data = data3Short,
							cluster = "justice",
							family = binomial(link = logit))
							
summary(model3VACBFullPooled$glm_res)
logLik(model3VACBFullPooled$glm_res)
BIC(model3VACBFullPooled$glm_res)
nobs(model3VACBFullPooled$glm_res)

models3 <- list(model2MeritsFullPooled$glm_res, model3MeritsFullPooled$glm_res, model2CVSGFullPooled$glm_res, model3CVSGFullPooled$glm_res, model2VACBFullPooled$glm_res, model3VACBFullPooled$glm_res)
modelsummary(models3, output = "latex", stars = TRUE)

#########################

###
# Control for first term new admin 
###

model5MeritsFullPooled <- glm.cluster(votesWithSG ~ policyChange
							+ noPolicyChange
							+ coeNoDetail
							+ coeCriminalCase
							+ coeInPrevious
							+ coePastTally
							+ realOutlierSignal
							+ firstTermNewAdmin
							+ sgJusticeDist
							+ sgPetitioner
							+ constitutionalCase
							+ civilRightsCase
							+ declareUncon
							+ CSI,
							data = data1,
							cluster = "justice",
							family = binomial(link = logit))
							
summary(model5MeritsFullPooled$glm_res)
logLik(model5MeritsFullPooled$glm_res)
BIC(model5MeritsFullPooled$glm_res)
nobs(model5MeritsFullPooled$glm_res)

###
# control for GVR dissent 
###

model6MeritsFullPooled <- glm.cluster(votesWithSG ~ policyChange
							+ noPolicyChange
							+ coeNoDetail
							+ coeCriminalCase
							+ coeInPrevious
							+ coePastTally
							+ realOutlierSignal
							+ dissentGVR
							+ sgJusticeDist
							+ sgPetitioner
							+ constitutionalCase
							+ civilRightsCase
							+ declareUncon
							+ CSI,
							data = data1,
							cluster = "justice",
							family = binomial(link = logit))
							
summary(model6MeritsFullPooled$glm_res)
logLik(model6MeritsFullPooled$glm_res)
BIC(model6MeritsFullPooled$glm_res)
nobs(model6MeritsFullPooled$glm_res)

###
# Control for first term new admin 
###

model5CVSGFullPooled <- glm.cluster(votesWithSG ~ policyChange
							+ noPolicyChange
							#+ coeNoDetail
							+ coeCriminalCase
							+ coeInPrevious
							+ coePastTally
							+ realOutlierSignal
							+ firstTermNewAdmin
							+ sgJusticeDist
							+ sgPetitioner
							+ constitutionalCase
							+ civilRightsCase
							+ declareUncon
							+ CSI,
							data = data2,
							cluster = "justice",
							family = binomial(link = logit))
							
summary(model5CVSGFullPooled$glm_res)
logLik(model5CVSGFullPooled$glm_res)
BIC(model5CVSGFullPooled$glm_res)
nobs(model5CVSGFullPooled$glm_res)

###
# control for GVR dissent 
###

model6CVSGFullPooled <- glm.cluster(votesWithSG ~ policyChange
							+ noPolicyChange
							#+ coeNoDetail
							+ coeCriminalCase
							+ coeInPrevious
							+ coePastTally
							+ realOutlierSignal
							+ dissentGVR
							+ sgJusticeDist
							+ sgPetitioner
							+ constitutionalCase
							+ civilRightsCase
							+ declareUncon
							+ CSI,
							data = data2,
							cluster = "justice",
							family = binomial(link = logit))
							
summary(model6CVSGFullPooled$glm_res)
logLik(model6CVSGFullPooled$glm_res)
BIC(model6CVSGFullPooled$glm_res)
nobs(model6CVSGFullPooled$glm_res)

###
# Control for first term new admin 
###

model5VACBFullPooled <- glm.cluster(votesWithSG ~ policyChange
							+ noPolicyChange
							+ coeNoDetail
							+ coeCriminalCase
							+ coeInPrevious
							+ coePastTally
							+ realOutlierSignal
							+ firstTermNewAdmin
							+ sgJusticeDist
							+ sgPetitioner
							+ constitutionalCase
							+ civilRightsCase
							+ declareUncon
							+ CSI,
							data = data3,
							cluster = "justice",
							family = binomial(link = logit))
							
summary(model5VACBFullPooled$glm_res)
logLik(model5VACBFullPooled$glm_res)
BIC(model5VACBFullPooled$glm_res)
nobs(model5VACBFullPooled$glm_res)

###
# control for GVR dissent 
###

model6VACBFullPooled <- glm.cluster(votesWithSG ~ policyChange
							+ noPolicyChange
							+ coeNoDetail
							+ coeCriminalCase
							+ coeInPrevious
							+ coePastTally
							+ realOutlierSignal
							+ dissentGVR
							+ sgJusticeDist
							+ sgPetitioner
							+ constitutionalCase
							+ civilRightsCase
							+ declareUncon
							+ CSI,
							data = data3,
							cluster = "justice",
							family = binomial(link = logit))
							
summary(model6VACBFullPooled$glm_res)
logLik(model6VACBFullPooled$glm_res)
BIC(model6VACBFullPooled$glm_res)
nobs(model6VACBFullPooled$glm_res)


models4 <- list(model5MeritsFullPooled$glm_res, model6MeritsFullPooled$glm_res, 
model5CVSGFullPooled$glm_res, model6CVSGFullPooled$glm_res, model5VACBFullPooled$glm_res,  model6VACBFullPooled$glm_res)
modelsummary(models4, output = "latex", stars = TRUE)

#################

###
# add fixed effects by administration 
###
# from BOWW 2020

model7MeritsFullPooled <- glmer(votesWithSG ~ policyChange
							+ noPolicyChange
							+ coeNoDetail
							+ coeCriminalCase
							+ coeInPrevious
							+ coePastTally
							+ politicization
							+ realOutlierSignal
							+ sgJusticeDist
							+ sgPetitioner
							+ constitutionalCase
							+ civilRightsCase
							+ declareUncon
							+ CSI
							+ (1 | presName),
							data = data1,
							family = binomial(link = logit))
							
summary(model7MeritsFullPooled)
logLik(model7MeritsFullPooled)
BIC(model7MeritsFullPooled)
nobs(model7MeritsFullPooled)

###
# fixed effects by president
###

model7CVSGFullPooled <- glmer(votesWithSG ~ policyChange
							+ noPolicyChange
							#+ coeNoDetail
							+ coeCriminalCase
							+ coeInPrevious
							+ coePastTally
							+ politicization
							+ realOutlierSignal
							+ sgJusticeDist
							+ sgPetitioner
							+ constitutionalCase
							+ civilRightsCase
							+ declareUncon
							+ CSI
							+ (1 | presName),
							data = data2,
							family = binomial(link = logit))
							
summary(model7CVSGFullPooled)
logLik(model7CVSGFullPooled)
BIC(model7CVSGFullPooled)
nobs(model7CVSGFullPooled)

###
# fixed effect for president
###

model7VACBFullPooled <- glmer(votesWithSG ~ policyChange
							+ noPolicyChange
							+ coeNoDetail
							+ coeCriminalCase
							+ coeInPrevious
							+ coePastTally
							+ politicization
							+ realOutlierSignal
							+ sgJusticeDist
							+ sgPetitioner
							+ constitutionalCase
							+ civilRightsCase
							+ declareUncon
							+ CSI
							+ (1 | presName),
							data = data3,
							family = binomial(link = logit))
							
summary(model7VACBFullPooled)
logLik(model7VACBFullPooled)
BIC(model7VACBFullPooled)
nobs(model7VACBFullPooled)

stargazer(model7MeritsFullPooled, model7CVSGFullPooled, model7VACBFullPooled)

###
# cluster standard errors by term
###

model8MeritsFullPooled <- glm.cluster(votesWithSG ~ policyChange
							+ noPolicyChange
							+ coeNoDetail
							+ coeCriminalCase
							+ coeInPrevious
							+ coePastTally
							+ politicization
							+ realOutlierSignal
							+ sgJusticeDist
							+ sgPetitioner
							+ constitutionalCase
							+ civilRightsCase
							+ declareUncon
							+ CSI,
							data = data1,
							cluster = "term",
							family = binomial(link = logit))
							
summary(model8MeritsFullPooled$glm_res)
logLik(model8MeritsFullPooled$glm_res)
BIC(model8MeritsFullPooled$glm_res)
nobs(model8MeritsFullPooled$glm_res)


###
# cluster standard errors by term
###

model8CVSGFullPooled <- glm.cluster(votesWithSG ~ policyChange
							+ noPolicyChange
							#+ coeNoDetail
							+ coeCriminalCase
							+ coeInPrevious
							+ coePastTally
							+ politicization
							+ realOutlierSignal
							+ sgJusticeDist
							+ sgPetitioner
							+ constitutionalCase
							+ civilRightsCase
							+ declareUncon
							+ CSI,
							data = data2,
							cluster = "term",
							family = binomial(link = logit))
							
summary(model8CVSGFullPooled$glm_res)
logLik(model8CVSGFullPooled$glm_res)
BIC(model8CVSGFullPooled$glm_res)
nobs(model8CVSGFullPooled$glm_res)


###
# cluster standard errors by term
###

model8VACBFullPooled <- glm.cluster(votesWithSG ~ policyChange
							+ noPolicyChange
							+ coeNoDetail
							+ coeCriminalCase
							+ coeInPrevious
							+ coePastTally
							+ politicization
							+ realOutlierSignal
							+ sgJusticeDist
							+ sgPetitioner
							+ constitutionalCase
							+ civilRightsCase
							+ declareUncon
							+ CSI,
							data = data3,
							cluster = "term",
							family = binomial(link = logit))
							
summary(model8VACBFullPooled$glm_res)
logLik(model8VACBFullPooled$glm_res)
BIC(model8VACBFullPooled$glm_res)
nobs(model8VACBFullPooled$glm_res)

models5 <- list(model8MeritsFullPooled$glm_res, model8CVSGFullPooled$glm_res, model8VACBFullPooled$glm_res)
modelsummary(models5, output = "latex", stars = TRUE)